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Abstract 

For the minimization of a nonlinear cost functional j under convex con¬ 
straints the relaxed projected gradient process 

ipk+i = <Pk + a k (PH(<Pk - A kVHj(<Pk)) ~ <Pk) 

as formulated e.g. in [ 12] is a well known method. The analysis is classically 
performed in a Hilbert space H. We generalize this method to functionals 
j which are differentiable in a Banach space. Thus it is possible to perform 
e.g. an L 2 gradient method if j is only differentiable in L°°. We show global 
convergence using Armijo backtracking in a k and allow the inner product and 
the scaling A*, to change in every iteration. As application we present a struc¬ 
tural topology optimization problem based on a phase held model, where the 
reduced cost functional j is differentiable in H 1 nL°°. The presented numerical 
results using the H 1 inner product and a pointwise chosen metric including 
second order information show the expected mesh independency in the iter¬ 
ation numbers. The latter yields an additional, drastic decrease in iteration 
numbers as well as in computation time. Moreover we present numerical re¬ 
sults using a BFGS update of the H 1 inner product for further optimization 
problems based on phase held models. 


Key words: projected gradient method, variable metric method, convex con¬ 
straints, shape and topology optimization, phase field approach. 
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1 Introduction 

Let j be a functional on a Hilbert space H with inner product (., ,)h and induced 
norm ||.||# and let <& a d ^ H be a non-empty, convex and closed subset. We consider 
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the optimization problem 


minj(^) subject to ip e <h a rf. 


( 1 ) 


If j is Frechet differentiable with respect to ||.||ff, the classical projected gradient 
method introduced in Hilbert space in [TSj and [23] can be applied, which moves in 
the direction of the negative //-gradient -V uj e //, which is characterized by the 
equality (V Hj(<^),v)H = {j'(<p)i v)h* h e H an d orthogonally projects the result 
back on to stay feasible, i.e. 


(fik+i = - X k y H j(^ k ))- 


( 2 ) 


To obtain global convergence A*, has to be chosen according to some step length rule, 
which results in a gradient path method, or one can perform a line search along the 
descent direction Vk = Pii(<Pk - XkVHj(<Pk)) ~<Pk- A typical application is H = L 2 (fi), 
see e.g. El 

In this paper we consider the case that j is differentiable with respect to a norm 
which is not induced by a inner product. Hence no //-gradient V jjj exists. However, 
in Section [2] we reformulate the method such that it is well defined under weaker 
conditions. We show global convergence when Armijo backtracking is applied along 
Vk and allow the inner product and the scaling Xk to change in every iteration. We 
call this generalization ‘variable metric projection’ type (VMPT) method. In Section 
[3] we study the applicability of the method to a structural topology optimization 
problem, namely the mean compliance minimization in linear elasticity based on a 
phase held model. Then the reduced cost functional is differentiable only in H 1 r\L°°. 
In the last section we show numerical results for this mean compliance problem. As 
expected choosing the H 1 metric leads to mesh independent iteration numbers in 
contrast to the L 2 metric. We also present the choice of a variable metric using 
second order information and the choice of a BFGS update of the H 1 metric. This 
reduces the iteration numbers to less than a hundreth. Moreover, we give additional 
numerical examples for the successful application of the VMPT method. These 
include a problem of compliant mechanism, drag minimization of the Stokes how 
and an inverse problem. 

2 Variable metric projection type (VMPT) method 

2.1 Generalization of the projected gradient method 

The orthogonal projection Pn{^Pk ~ A^V Hj{^k)) employed in (J2]) is the unique solu¬ 
tion of 
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which is equivalent to the problem 

min \ \y - (p k \\ 2 H + A k Dj(<p k ,y- <p k ), (3) 

y^ad 2 


since (V H j(p k ),y ~ Vk)H = j'{yk){y ~ Pk) = Dj((p k ,y - p k ) where the last denotes 
the directional derivative of j at ip k in direction y - Lp k . If e.g. Dj(tp k ,y) is linear 
and continuous with respect to y € H the cost functional of (J3]) is strictly convex, 
continuous and coercive in H, and hence ([3]) has a unique solution (p k nm- In the 
formulation ([3]) the existence of the gradient V hJ is not required. Even Gateaux 
differentiability can be omitted. 

In the following we formulate an extension of the projected gradient method where 
- Ais replaced by the solution (p k of (|3]). 

First we drop the requirement of a gradient as mentioned above. We assume that 
the admissible set Q a d is a subset of an intersection of Banach spaces InD, where 
X and D have certain properties (see (Al) ), which are e.g. fulfilled for X = i/ 1 (G) 


or X = L 2 {yi) and B = L°°{yi). Furthermore assume that j is continuously Frechet 
differentiable on d> aa f with respect to the norm ||.||xnro : = l-lx + 11 -Id- The Frechet 
derivative of j at <p is denoted by j'(p) e (X n D)* and we write (.,.) for the dual 
paring in the space X n B. Moreover, we use C as a positive universal constant 
throughout the paper. 

Secondly, we also allow the norm ||.||# in (J3| to change in every iteration. Therefore, 
we consider a sequence {a k } k >o of symmetric positive definite bilinear forms inducing 
norms ||.|| afc onXnD . This approach falls into the class of variable metric methods 
and includes the choice of Newton and Quasi-Newton based search directions (see for 
example mm and [T9] for the unconstrained case). In j2j these methods are called 
scaled gradient projection methods and in the case of a k = j"(<Pk) also constrained 
Newton’s method. In finite dimension a k is given by a k (p,v ) : =p T B k v where B k can 
be the Hessian at ip k or an approximation of it. 

Hence, in each step of the VMPT method the projection type subproblem 


min h\y-<p k \\l k + ^ k {j'(Vk),y-<Pk) (4) 

y^ad 2 k 

with some scaling parameter A& > 0 has to be solved. Problem (jd]) is formally 
equivalent to the projection P ak (tp k - A k Va k j{<£k))- However, j is not necessarily 
differentiable with respect to ||.|| afc and XnD endowed with a k (.,.) is only a pre- 
Hilbert space. Hence V ak j{p k ) does not need to exist. For globalization of the 
method we perform a line search based on the widely used Armijo back tracking, 
which results in Algorithm |2.1| In the next section it is shown that the algorithm is 
well defined under certain assumptions and in particular that a unique solution (p k 
of (Jd]) exists, together with the proof of convergence. We denote the solution of Q 
also by V k ((p k ) due to the connection to a projection. 


3 





Algorithm 2.1 (VMPT method). 

1: Choose 0 < /3 < 1, 0<cr<l and tpo e <$> a( i. 

2: k:= 0 

3: while k < k max do 
4: Choose A fc and a k . 

5: Calculate the minimum Tp k = V k (vk) of the subproblem Q. 

6: Set the search direction v k Vk - cp k 

7: */||v*||x£ tol then 

8: return 

9: end if 

10: Determine the step length Oik '■= /3 mk with minimal rrik e N 0 such that 

j(Vk + Oi k v k ) <j(Vk) + a k o(j'(ip k ),v k ). 

11: Update ip k +1 := Vk + OL k v k 

12: k := k + 1 

13: end while 


The stopping criterion ||ufc|x < tol is motivated by the fact that ip k is a stationary 
point of j if and only if v k = 0 and v k -*■ 0 in X, cf. Corollary |2.6 and Theorem T2 


We would like to mention, that this algorithm is not a line search along the gradient 
path , which is widely used (e.g. in (2 []J. [El El US, US ED, EH E5]) and which 
requires to solve a projection type subproblem like (|2]) in each line search iteration. 
This can be unwanted if calculating the projection is expensive compared to the 
evaluation of j. To avoid this we perform a line search along the descent direction 
v k , which is suggested e.g. in finite dimension or in Hilbert spaces in PIESEl] and 
is also used in [[13]. To include the idea of the gradient path approach, we imbed 



the possibility to vary the scaling factor {\ k } k >o for the formal gradient in 
each iteration. The parameter X k can be put into a k by dividing the cost in (|4j) 

Afc. However, we treat it as a separate parameter since this reflects the case where 
a k is fixed for all iterations. Note that under the assumptions used in this paper a 
line search along the gradient path is not possible since not even the existence of a 
positive step length can be shown, cf. Remark 2.8] 


Moreover, there is a clear connection to sequential quadratic programming, consid¬ 
ering that VkisPk) is the solution of the quadratic approximation of min ¥ , e $ ad j(<p) 
with 

min j(<p k ) + {f(<p k ),y- Vk) + J a k (y -<p k ,y- Vk)- 

y^ad 2 

However, the global convergence result is analysed by means of projected gradient 
theory. 
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2.2 Global convergence result 

We perform the analysis of the method with respect to two norms in the spaces X 
and ID, which we assume to have the following properties: 

(Al) X is a reflexive Banach space. ID is isometrically isomorphic to IB*, where B 
is a separable Banach space. Moreover, for any sequence {<p;} in X n B with 
ifi -> <p weakly in X and ipi -*■ (p weakly-* in ID, it holds ip - (p. 

We identify ID and B* and say that a sequence converges weakly-* in ID if it converges 
weakly-* in B*. The separability of B is used to get weak-* sequential compactness 
in ID. We would like to mention that the results hold also if ID is a reflexive Banach 
space, in particular if ID is an Hilbert space. In this case weak-* convergence has to 
be replaced by weak convergence throughout the paper. However, in the application 
we are interested in ID = 

In case of the Sobolev space X = W k ' p {VL) and ID = where H is a bounded 

domain, k > 0, 1 < p < oo and 1 < q < oo the above assumption is fulfilled. 

In addition to the above conditions on X and ID let the following assumptions hold 
for the problem (JT|: 

(A2) ‘had £ X n ID is convex, closed in X and non-empty. 

(A3) <h a d is bounded in ID. 

(A4) j(<p) > -C > -oo for some C > 0 and all <p e <h ad . 

(A5) j is continuously differentiable in a neighbourhood of < h a( i <= X n ID. 

(A6) For each ip e <h a d and for each sequence {pi] [InD with ipi -*■ 0 weakly in 
X and weakly-* in ID it holds {j'(<p),(pi) -*■ 0 as i -> oo. 

Moreover, we request for the parameters a k and A*, of the algorithm that: 

(A7) {a k } is a sequence of symmetric positive definite bilinear forms on X n ID. 

(A8) It exists c\ > 0 such that ci|p|^ < \\p\\l k for all p e X n B and k e N 0 . 

(A9) For all k 6 No it exists C 2 (At) such that < C2||p||xnro ^ or all p 6 X n D. 

(A10) For all k 6 N 0 , p e and for each sequence {i/j} £ <h a( i where there exists 
some i/eXnID with y t -*■ y weakly in X and weakly-* in B it holds a k (p,yi) 
a k (p,y) as z ^ oo. 

(All) For each subsequence {<p ki }i of the iterates given by Algorithm [2H] converg¬ 
ing in X n ID, the corresponding subsequence {a ki }i has the property that 
dkiiPiiVi) 0 for any sequences {pi}, {yi] £ X n ID with pi -> 0 strongly in X 
and weakly-* in ID and {z/j} converging in X n ID. 

(A12) It holds 0 < A m i n < A k < A„ for all k € N 0 . 
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(Al) (A12) are assumed throughout this paper if not mentioned otherwise. 


Assumption (All) reflects the possibility of a point based choice of a,k, e.g. de¬ 
pendent on the Hessian D 2 j(ipk) or on an approximation of the Hessian. Note that 

< c 2 |bl 2 


(A9) (All) is weaker than the assumption ||p|| 


2 

a k 


IIX- 


Ofc is given, which only fulhlls these weaker assumptions. Also (A8) is weaker than 


In (21) an example of 


C l|l M llxnjp_ 

Section 12.31 


< ||u||^ . The main result of the paper is the following, which is proved in 


Theorem 2.2. Let {</?&} £ &ad be the sequence generated by the VMPT method 
(Algorithm \2.1\ ) with tol = 0 and let the assumptions (Al) (A 12) hold, then: 

1. lim k ^ooj(Pk) exists. 

2. Every accumulation point of {pk) in X n D is a stationary point of j. 

3. For all subsequences with <p fci -»■<£> inXnD where p is stationary, the subse¬ 
quence {vk, },, converges strongly in X to zero. 

4- If additionally j e C 1,7 ^^) with respect to ||.||xnro for some 0 < 7 < 1 then the 
whole sequence {vk}k converges to zero in X. 

In the classical Hilbert space setting, i.e. D = X = H for some Hilbert space H , the 


assumption (A3) can be dropped. Also assumption (A6) is trivial because of (A5) 


Moreover, assumptions KArn(Aii)| are fulfilled for the choice ak(p,v) = ( p,Akv)n 
where Ak e C(H) is a self-adjoint linear operator with m||p||# < ( p,Akp)h ^ 
and M > m > 0 independent of k. This is e.g. assumed in the local convergence 
theory in dam! and in finite dimension for global convergence in [21 [23]. For the 
special choice ak(p,v) = (p,v)n, global convergence is shown in [Tj5] and for the case 
of a line search along the gradient path in 0. Result 4. of Theorem |2.2| is shown 
in ra in case of a line search along the gradient path under the same assumption 
j € C 1,7 . Thus the presented method is a generalization of the classical method in 
Hilbert space. 

We would also like to mention the following: 

Remark 2.3. If there exists C > 0 such that [|p||b ^ Cj|p||x for all p £ 3L r D, 


assumption (A3) can be omitted . 


IfJLisa Hilbert space, the choice ak(u,v ) = ( u,v)h fulfills all assumptions (Al)- 

TAm 


2.3 Analysis and proof of the convergence result of the VMPT 
method 

We first show the existence and uniqueness of Tp k = Vk(p>k) based on the direct 
method in the calculus of variations using the following Lemma and assumptions 
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(Ml EMI and [(A5j]|(A10)1 Note that the standard proof cannot be applied, 


since a k is indeed X-coercive, but a*, and (j'(p k ),-) are not X-continuous. Another 
difficulty is that X n D is not necessarily reflexive. 


Lemma 2.4. Let {p k } £ < h a d with p k -*■ p weakly in X for some p e <F a d- Then p k -> p 
weakly-* in D. 


Proof. Since <F a d is bounded in D and the closed unit ball of D is weakly-* sequen¬ 
tially compact due to the separability of B, we can extract from any subsequence 
of {p k } £ $ad another subsequence {p ki } with p k . -*■ p weakly-* in D for some peD. 
Due to the required unique limit in X and D we have p = p. Since for any sub¬ 
sequence we find a subsequence converging to the same p, we have that the whole 
sequence converges to p. □ 


Theorem 2.5. For any k e No and p e <P a d, the problem 

min ^\\y-(p\\l h + X k {j'(p),y-p) (5) 

admits a unique solution p := V k (p), which is given by the unique solution of the 
variational inequality 


a k {p-p,iq-p) + \ k {j'(p),ri-p) >0 Mr) e <f> ad . 

Proof. Let k e N 0 and p € <f> ad arbitrary. Problem (|5]) is equivalent to 

min g k (y) := \a k (y,y) + {b k ,y) 
ye&ad 


( 6 ) 

(7) 


where {b k ,y) := X k (j'(p),y) - a k (p,y) and b k e (XnD)‘ due to (A5) and (A9) By 


(A3) and |(A8) we get for any y e <F a d with some generic C > 0 

I b k | (XnD)* (||?/||x + IMIb) k ~C. 


9k(y) > y||j/||x 


( 8 ) 


<c 


Thus g k is X-coercive and bounded from below on <J> a d. Hence we can choose an 
inffinizing sequence pi € <h arf , such that g k (pi ) °°> inf y6 $ od g k (y). From the estimate 


we conclude that {pi}i is bounded in X. Therefore, we can extract a subsequence 
(still denoted by pf) which converges weakly in X to some (p e X. Since 'Fad is convex 
and closed in X, it is also weakly closed in X and thus p € <F a d- By Lemma 2 A we also 
get pi ^ p weakly-* in ID. Finally we show g k (p) = inf ye $ od g k (y). Using|(A6)j|(A8) 


and (A10) one can show that liminfj a k (p i: pf) > a k (p,p) and lim, (b k ,pi) = (b k ,p), 
thus liminfj g k {pf) ^gk(p)- We conclude 


inf gk(y) < gk(p) < lim.infpjfc(^i) = inf g k (y), 

y^ad * y^ad 
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which shows the existence of a minimizer of (JT]) . Using (A 8 ) the uniqueness follows 
from strict convexity of g k . 

Due to KAgfl and |(A9)j we have that gk is differentiable in X n ID, where its direc¬ 
tional derivative at p in direction rj - p for arbitrary rj e Q ad is given by 


{a'M’V-v) = a k (p-p,rj-p) + X k {j'(p),g-p) . 

Since the problem ([5]) is convex, it is equivalent to the first order optimality condi¬ 
tion, which is given by the variational inequality (Jb]), see [55]. □ 

We see that p e & ad is a stationary point of j, i.e. {j'(p),ij - (p) > 0 V 77 e & ad , if 
and only if p - p is the solution of (J 6 ]) , i.e. the fixed point equation p = V k (p) is 
fulfilled. This leads to the classical view of the method as a fixed point iteration 
<Pk+i = Vk(<Pk) in the case that V k is independent of k and a k - 1 is chosen. 

Corollary 2.6. If there exists some k € N 0 with V k (p) ~ T then ip is a stationary 
point of j. On the other hand, if p € Q ad is a stationary point of j then the fix point 
equation V k fp) = p holds for all k e No- In particular, an iterate p k of the algorithm 
is a stationary point of j if and only if v k = Vk(pk) - p k = 0. 

The variational inequality (| 6 ]) tested with rj - p € & ad together with 
yields that V k (p) - <p is a descent direction for j: 

Lemma 2.7. Let k e No, p e <& ad and v Vk{p) - p- Then it holds 

( ?(p),v) <-T^lk|||. (9) 

^max 

□ 


(A 8 ) and (A12) 


Note that ([9]) does not hold in the X n D-norm. 

Due to ( j'(p),v) < 0 for v + 0 the step length selection by the Armijo rule (see step 
10 in Algorithm |2.l[) is well defined, which can be shown as in |2J. 


Remark 2.8. For the existence of a step length and for the global convergence proof 
we exploit that the path a p k + av k is continuous inln D. Thus, also the mapping 
a i-» j(pk+OiV k ) is continuous. On the other hand, this does not hold for the gradient 
path. Backtracking along the gradient path or projection arc means that a k is set to 
1, whereas \ k = /3 mk is chosen with m k € Nq minimal such that the Armijo condition 


iivki^k)) < j{pk) + o {j\p k ), p k (\ k ) - p k ) 


is satisfied, see for instance jMf. By the notation p k (\ k ) we emphasize that the 
solution of the subproblem Q depends on \ k . However, with the above assumptions 
it cannot be shown that there exists such a \ k . The reason is that due to 17^1 
the gradient path A p k (X) is continuous with respect to the IL-norm, whereas j is 
due to (A5) only differentiable with respect to the InD -norm. Thus, j along the 
gradient path, i.e. the mapping X j(p k ( A)), may be discontinuous. 















To prove statement 2. of Theorem 2.2 we use, as in [2] for finite dimensions, that 
Vk is gradient related. This is weaker than the common angle condition. Therefor 
we need the following two lemmata: 

Lemma 2.9. For {(fik}k - &ad with ipk ^ (f mXn B and {pk}k £XnD with pk -*■ p 
weakly in X and weakly-* in B for some <p,pz XnB it holds {j'((fk),Pk) (j'i^iP)- 

Proof. We use EMI and |(A6)| and obtair 

I {j'(Vk),Pk) - {j'(p),p) I < I (j'ivk) - j'{v),Pk) I + I (. j'(<p),Pk-P) I ^ 

^ \\j'(<Pk) -/(vOlkxnD)* \\pkhnD + \{3'(p),Pk-p)\^0. □ 

'-V--V-' "-V-' 

->■0 <c ^0 


The preceding lemma is also needed in the proof of Theorem 2.2 


Lemma 2.10. Let for a sequence {pi]i £ ‘had hold ipi -»■ tp in X n B for some 
^eXnB. Then there exists C > 0 such that \\V k (Pi) ||xnro < C for all i,k e N 0 . 

Proof. Lemma 2.7 yields together with |( A3) and EM the estimate 


Cl 


■\\^k(Pi) - Ti\\x £-{j , (<Pi),'P k (<Pi)-<Pi) 


- \\j' (Pi)\\ (X.r\lD)* (\\'Pk( ( Pi) ~ 7h||x + || 'Pk(Ti) ~ Thill)) 

< C(\\Vk(<Pi) - <^i||x + 1), 


thus | Vk(<fi) ~ <£>t||x ^ C and hence \\Vk(pi)\\ji < C. Due to (A3) we finally get 
W'Pk(Ti )||xnD ^ C independent of i and k. □ 


Lemma 2.11. Let {pk) be the sequence generated by Algorithm 2.1. then {i’k}k is 
gradient related, i.e.: for any subsequence {pk,}t which converges m XnB to a 
nonstationary point p e Q a d of j, the corresponding subsequence of search directions 
{vki}i Is bounded in XnB and limsupj {j'(Pki)i v ki) < 0 is satisfied. Moreover, it 
holds lim inf j |k^ ||x > 0. 


Proof. Let p ki ->^inXnB, where p is nonstationary. Lemma 2.10 provides that 
{v^i is bounded in XnB. With Q, the statement limsup, (j'(p k ,f), Vk.,) < 0 follows 
from lim inf', ||nfcj|x = C > 0, which we show by contradiction. 

Assume lim inf* fvfcjx = 0, thus there is a subsequence again denoted by {v ki }i such 
that Vk, -*■ 0 in X. Using (J(3|) for pk ■= Pk(Tk), the positive definiteness of ak and 
IAI2)I it follows for all rj e <h a d 


{j'(.Tk),V -p k )> t r( a k(v k , v k ) + a k (v k , p k - v k -rj)) 

''K 

- -\Z-\ a k( v k,Pk-v k -v)\ ■ 


( 10 ) 
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Moreover, p k . = v k . + p kt -*■ <p in X and also weakly-* in D according to Lemma 


2.4. From Lemma 2.9 we get (j'(p ki ),rj ~ pki) -*■ (j'(<p), V ~ <p)- From (All) we get 


a k t (<Pki ~ <Pki,<Pki - V) -*• 0 and we derive from (|T0j) that 

0 Vr7€$ ad , 

which shows that p is stationary, which is a contradiction. 


□ 


Proof of Theorem 2.2 


Because of Corollary |2.6| we can assume v k y 0 and a k > 0 for all k. 

1.) From the Armijo rule and since v k is a descent direction we get 

j(<Pk+i) ~ j(<Pk) ^ OL k a{j'(ip k ),v k ) < 0, 


( 11 ) 


thus j(p k ) is monotonically decreasing. Since j is bounded from below we get 
convergence j{p k ) j* f° r some j* e R, which proves 1. 

2.) The proof is similar to [2] in finite dimension by contradiction. Let p be an 
accumulation point, with a convergent subsequence p ki T i n XnB. The continuity 
of j on <F a(i yields then j* = j{p) and (11) leads to a k (j'(p k ),v k ) 0. Assuming 


now that p is nonstationary we have \{j'(pk,,), v k,)\ > C > 0, since {v k } is gradient 
related by Lemma 2.11, and thus a ki -> 0. So there exists some j e N such that 
aikjP < 1 for all i > i, and thus akjft does not fulfill the Armijo rule due to the 
minimality of m*.. Applying the mean value theorem to the left hand side, we have 
for some nonnegative ay. < and all i > i that 

-J ■ (f i'Pki + , v k f) = j [p K + ^ v ki J - j(p ki ) > ^a {j'(p ki ),v ki ) (12) 


holds. Since, by Lemma 2.11 

Vki+Oi ki v ki 


{v ki }i is bounded in X n D and a ki ->• 0, we have that 


(^inXnD. Also p k . = p k .+v ki is uniformly bounded inXnD and thus 
there exists a subsequence, again denoted by {p ki }, which converges to some y e <F a(i 
weakly in X and weakly-* in D. Hence we have that v ki = p k . - p k . -»• v y - p 
weakly in X and weakly-* in D. According to Lemma 2.9 we can take the limit of 


both sides of the inequality (12), which leads to ( j' (p) ,v) > a (j’ (p) ,v ), and a < 1 
yields {j' (p) ,v) > 0 . This contradicts (j' (p), v) = limsup i(j'(p ki ),v ki ) < 0, which 


is a consequence of Lemma 2.11 


3.) By proving that out of any subsequence of ( j r (p ki )i v ki) we can extract another 
subsequence, which converges to 0, we can conlude that {j'(p ki ),v ki ) -* 0 which 
yields fryjlx ->■ 0 by (g. With Lemma |2.10[ we get by the same arguments as in 
2. that v ki -»• y - p weakly in X and weakly-* in D for a subsequence and for some 


y e $ ad , thus ( j'{pki),Vk t ) -*■ (j'(p),y - p) due to Lemma 2.9. Since v ki are descent 


directions for j at p ki and p is stationary we have ( j'{p),y - p) - 0. 


10 
















4.) As in 3.) we prove by a subsequence argument that {j r ((fk)i v k) 0- For an 


arbitrary subsequence, which we also denote by index k, ( 11 ) yields a k {j'(pk),Vk) 

0. If a.k > c > 0 for all k, the assertion follows immediately. Otherwise there exists 
a subsequence (again denoted by index k) such that j3 > a k ^ 0 and thus the step 
length otk/P does not fulfill the Armijo condition. Since j' is Holder continuous with 
exponent 7 and modulus L we obtain 




{j'(<Pk),Vk) <j(vk + jrv k )-j((p k ) = J o £i( ( Pi< + t T Vk ') dt 
- U'iP^iVk) + (^-) IkfcllxnID- 

It holds |ufc|iD ^ C due to (A3) and employing (|9]) we obtain 

0 < (a - 1) <j'fe),»*> < Cj^(f rdl^llt* 7 + 1) < Cal (I (/(wO.**) I ¥ + 1). 

We get Xk ’■= | {j'(Tk),v k ) | ->• 0. Otherwise there exists a subsequence still denoted by 

- 1+7 

{xk} with Xk -*■ c > 0. Rearranging the last inequality gives 1 < Cof(x k 2 ->■ 0, 


which is a contradiction. 


□ 


Remark 2.12. Statements 1. and 2. of Theorem 2.2 require only that lp k e <3> a d 
is chosen such that the search directions v k = Tp k - ip k are gradient related descent 
directions, as can be seen in the proof above. Hence Tp k does not have to be Vkfpk) 


in Algorithm 2.1 In this case assumption (A3) is also not required . 


3 An application in structural topology optimiza¬ 
tion based on a phase field model 

In this section we give an example of an optimization problem described in [4], which 
is not differentiable in a Hilbert space, so the classical projected gradient method 
cannot be applied, but the assumptions for the VMPT method are fulfilled. 

We consider the problem of distributing N materials, each with different elastic 
properties and fixed volume fraction, within a design domain H £ R d , d e N, such 
that the mean compliance f r g u is minimal under the external force g acting on 
T„ 7 The displacement field u : -»• is given as the solution of the equations 


of linear elasticity (14). To obtain a well posed problem a perimeter penalization 
is typically used. Using phase fields in topology optimization was introduced by 
Bourdin and Chambolle [8] . Here, the N materials are described by a vector valued 
phase held ip : H -> R^ with p > 0 and ffi Ti = 1, which is able to handle topological 
changes implicitly. The ?’th material is characterized by {qp i = 1} and the different 
materials are separated by a thin interface, whose thickness is controlled by the 
phase held parameter e > 0. In the phase held setting the perimeter is approximated 
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by the Ginzburg Landau energy. In [5] it is shown that the given problem for N = 2 
converges as e -*■ 0 in the sense of T-convergence. For further details about the 
model we refer the reader to j4]J. The resulting optimal control problem reads with 

E (<p) := /o{|lv ^| 2 + - £ Mv)} 



min J(tp,u)-.= / g-u + ^E{ip) 

JVg 

(13) 


V e u e Hp ■— {H\Q) d | £| rD = 0} 


subject to 

f C(cp)£(u) :£(£)= f g.£ VZeHb 

Jn Jr g 

(14) 


r N 

f <p = ™, v >°, E^ e1 > 

Jn 

(15) 


where 7 > 0 is a weighting factor, <p := p f n <p, ip 0 '■ R /V -»• R is the smooth 
part of the potential forcing the values of ip to the standard basis e* e R v , and 
A : B := T,ij =1 AijBij for A, B e R dxd . The materials are fixed on the Dirichlet 
domain T £> n Q£l. The tensor valued mapping C : H N -*■ R dxd <£> (R dxd )* is a 
suitable interpolation of the stiffness tensors C(e{) of the different materials and 
£(u) \{yu + \/u T ) is the linearized strain tensor. The prescribed volume fraction 
of the zth material is given by trq. For examples of the functions -0 O an d C we refer to 
GUI- Existence of a minimizer of the problem (b3J) as well as the unique solvability 
of the state equation (14) is shown in [4] under the following assumptions, which we 
claim also in this paper. 


(AP) 17 n is a bounded Lipschitz domain; r D ,r g n QQ with T D r\T g = 0 and 
l-L d ~ l (T D ) > 0. Moreover, g e L 2 (Y g ) d and -00 e G 1 ’ 1 (R Ar ) as well as m > 
0, EiliUp = 1. For the stiffness tensor we assume C = (Cijki) d g kl=1 with 
Cijki e C 1 , 1 (R Ar ) and Cijki = CjM = Ckiij and that there exist ao,a±,C > 0, s.t. 
a 0 |A| 2 < C(cp)A : A < a 1 |A| 2 as well as |C»| < C holds for all symmetric 
matrices A e R dxd and for all <p e R w . 


The state u can be eliminated using the control-to-state operator S', resulting in 
the reduced cost functional j(<p) '■= J\<p , S (</?)). In |4] it is also shown that j : 
H 1 (£l) N n L°°(Q) N -> R is everywhere Frechet differentiable with derivative 

j'(<p)v= 7 AeVyj: W + -^(^M - f C'(ip)v£(u):£(u) (16) 

Jn e Jn 

for all tp,v e H 1 ( fl) N n L°°( £l) N , where u = S(<p ) and S ■ L°°( Q) N -»• H 1 ( Vt) d 
is Frechet differentiable. By the techniques in [4j one can also show that S' is 
continuous. 

In [H EJ the problem is solved numerically by a pseudo time stepping method with 
fixed time step, which results from an L 2 -gradient flow approach. An H~ l gradient 
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flow approach is also considered in [5J. The drawbacks of these methods are that 
no convergence results to a stationary point exist, and hence also no appropriate 
stopping criteria are known. In addition, typically the methods are very slow, i.e. 
many time steps are needed until the changes in the solution ip or in j are small. 
Here we apply the VMPT method, which does not have these drawbacks and which 
can additionally incorporate second order information. 

Since n L°°(Q) N is not a Hilbert space the classical projected gradient 

method cannot be applied. In the following we show that problem (13) fulfills the 
assumptions on the VMPT method. Amongst others we use the inner product 
a k(f,d) - Jq V/ : S7g. To guarantee positive definiteness of this a*, we first have 
to translate the problem by a constant to gain f } tp = 0 , which allows us to apply 
a Poincare inequality. Therefor we perform a change of coordinates in the form 
Cp - p - m and get the following problem for the transformed coordinates. 


mmj(p):= / g ■ S(p + m) + r yE(p + m) 

Jv g 

V s := Jv> <5 Jy 1 (f2) A ' f a V = 0, <p > -m, 2>- s oj. 


(17) 


On the transformed problem (17) we apply the VMPT method in the spaces 


X:= \p^H l {FL) 


N 


ih =0 }’ 


D := L°°(VL) 


N 


The space of mean value free functions X becomes a Hilbert space with the inner 
product (/,g)x : = (V/, Vg)m an d ||.||x is equivalent to the H 1 -norm [T|. 


Theorem 3.1. The reduced cost functional j : X n D 
differentiable and j' is Lipschitz continuous on <& a d- 


R is continuously Frechet 


Proof. The Frechet differentiability of j onXnD is shown in jlj. Let 77 , p i e Xn D 
and Ui = S(p j ), i = 1,2. Then with (16), ipo e C 1 , 1 (R JV ), Cijki e C' 1 , 1 (R JV ) and 
\C\p)\ < C \/<p e E w we get 


^ MWi - V2II m Mm + Cj\\tp 1 - p 2 \\ L 2 |t7|| L 2 

+ l/o( c '( m + <Pi) - C 'M + P 2 ))M£( u i) -SM)\ 

+ \L C '( xn + F 2 )M^( u l- u 2 ) :f(«l)| 

+ \ L C ' ( m + P2 )M£( u 2) ■S(Ui-U 2 )\ 

< C\\p l - p 2 \\ H i Mm 

+ |(C'(m + pf) - C\m + <p 2 )ML°°M\\m + 

+ CII^IU 00 1Wi - u 2 Iff 1 (IttiI// 1 + ||'u 2 |rr 1 ) 

< CMmrL°°{\\Pi - <p 2 ||m + Wi ~ 

+ ||iti - W2||r/ i (||'^i||H 1 + I u 2 1 m )} (18) 
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To show the continuity of j', let p n , e X n D for n e N with p n -»• p in X n D. 
Using (18) yields 


\\j'(Pn) ~ ?(<P)\hwnL~)* 

< C(\lp n - + |«n||tfi) + \\u n -u\\ H i(\\u n \\Hi + IttllHi)), 

where u n = S(p n ) and u = S(p). From the continuity of S we get that Iwnllrr 1 is 
bounded and that \u n - u||#i -»■ 0 as n -»■ oo. This implies 

II j'(Pn) -?(<P)hwnL~)* -0 


and thus j e C' 1 (X n D). 

For the Lipschitz continuity of j' we employ estimate (18) with p t e i = 1,2. 
Since <& a d is bounded in L°°, we get that S is Lipschitz continuous on <& a d and that 
\\S(v)\\m < C, independent of p e <3> a d, see [3j. This yields 


\j\Pl) -/(¥ 5 2 )ll(^ 1 nL“)* < C|| V 1 - V’all/finL-, 
which proofs the Lipschitz continuity of j' in <3> a d. 


□ 


Corollary 3.2. The spaces X and D, together with j and <f> a d given in ( [Tt| ) fulfill 
the assumptions MlWAeJ of the VMPT method. 


Proof. Given the choices for X and D (Al) is fulfilled. For p e & a d we have 

-1 < -m <^><l-m<l Vp e <3> ad 

almost everywhere in fh Thus it holds (A3) | and Q a d £ X n D. Moreover, 0 € $ QC j, 
<h a( i is convex, and since <1?^ is closed in L 2 (fii) N , it is also closed in X L 2 (Q) N . 

Thus (A2) holds. 


Assumption (A4) is shown in [1J and Theorem 3.1 provides (A5) 

Given 

if(p), Pi) = fjrt^P : + (^Mp + ™) - VC( P + m)£(u) : £(u)) ■ Pi} 


the first term converges to 0 if p { -*■ 0 weakly in H l . With (AP) and u e Hfi we 
have that o(p + m) - S/C(p + m)£(w) : S(u) e L 1 (h2) 7V . Hence the remaining 
term converges to 0 if p { -»• 0 weakly-* in L°°, which proves that (A6) is fulhlled. □ 


Possible choices of the inner product a*, for the VMPT method are the inner product 
on X, i.e. 


aifip^y) = (P,y)x = 



Vp: Vp 


(19) 
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and the scaled version a k (p,y) = 7 e(p,y)x- Both fulfill the assumptions (A7)- 


(All) We also give an example of a pointwise choice of an inner product, which 


includes second order information. Since this choice is not continuous in X, it is not 
obvious that it fulfils the assumptions. To motivate the choice of this inner product 
we look at the second order derivative of j, which is formally given by 


3 


i"( l fik)[p,y] = fjt7^P-Vy-2{C\m + ip k ){y)£{S\ip k )p):£{u k ))+ 

+ -V 2 ^ 0 (m + (p k )p ■ y - C"(m + cp k )[p,y]£(u k ) ■£(u k )}. 


In [4] it is shown that z p := S'(tp k )p e Hf is the unique weak solution of the linearized 
state equation 


f C(m + ^ k )£(z p ) :£(ri) = - [ C'(m + <p k )p£(u k ) : £{rj) \JpeH 
Jn Jn 


( 20 ) 


and that ||;2p||#i < C'||p||£°° holds. Since the first two terms in j" define an inner 
product (see proof of Theorem 3.3), we use 

ak(p,y)=1f£(p,y)x-2 / C'(m + cp k )(y)£(z p ):£(u k ) 

Jn 


( 21 ) 


as an approximation of j"(<p k ). Testing equation (20) for z y = S'((p k )y with z p we 
can equivalently write 




(p,y) = J£(P,y)x+2 J C(m+ <p k )£(z p ) ■■S(z t ). 


( 22 ) 


We would like to mention that the C 2 -regularity of j is not necessary for this defi¬ 
nition of a k . 


Theorem 3.3. The bilinear form a k given in (21) fulfills the assumptions (A7)~ 

Jam 1 


Proof. Due to (AP) and (22) we have 

ak(p,p) * T-bill- 

Thus, (A7) and (A8) is fulfilled. Furthermore, |( A9j holds due to 

a>k(p,y) + C\\z p \\ H i\z y \\ H i 

^ T e \\p b 1 WvWh 1 + C\\p\\l°° ||2/||l“ ^ C'lblUnroblUnD- 


(A10) is proved as in Corollary 3.2 


Finally we prove (All) For y k 0 and p k -*■ p in X we have (y k ,p k )x -»• 0 
for k -> 00. With <p k -»• <p, p k -> p in D = L°°(Q) N and S ■ L°°(Q) N -> if 1 (f]) JV 
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continuously Frechet differentiable, we have u k = S(tp k ) -*■ S(<p) =: u in H l D and 
z Pk = S'(<p k )p k -> S'(ip)p =: z p in Hjj. In particular, the sequences are bounded in 
the corresponding norms, including ||t/ fc ||£°° < C if y k -> y weakly-* in L°°. Using the 
Lipschitz continuity and boundedness of C' and VC(m + ip)£(z p ) ■ £(u) e L 1 (U) JV 
we have 


\LC'(m + cp k )y k £(z Pk ) : £(u k )\ 

^ \fn( C '( m + c Pk)-C'(m + cp))y k £(z Pk ) :£(u k )\ 

+ \L c '( m + ( p)yk£( z P k - z p) : £(uk)\ 

+ I L c '( m + v)Vk£ Op) -£(u k -u) | + | / n C \m + ip)y k £(z p ) : £{u)\ 

^ L\\<Pk _ vh-lVkh- \\ z Pk Wh 1 lOfclOi 
+ ||C"(m + lp)\l°° \\y k \\L°°\\z Pk -z p \\m\\ u k\\m 
+ |C'(m + ip) | L oo | y k | L oo | z p | H i | u k -u | H i 
+ l/n(VC'(m + 9?)£(z p ) :£{u))-y k \ ^ 0, 


which gives (All) 


□ 


Hence with 0 < A min ^ A k < A max-, all assumptions of Theorem |2.2| are fulhlled and 
we get global convergence in the space H l (£l) N n L°°(£i) N . 


4 Numerical results 


We discretize the structural topology optimization problem (13)-(15) using standard 


piecewise linear finite elements for the control Lp and the state variable u. The pro¬ 
jection type subproblem Q is solved by a primal dual active set (PDAS) method 
similar to the method described in |7|. Many numerical examples for this problem 
can be found in M, e.g. for cantilever beams with up to three materials in two 
or three space dimensions and for an optimal material distribution within an airfoil. 
In |3] the choice of the potential -0 as an obstacle potential and the choice of the 
tensor interpolation C is discussed. Also the inner products (., -)x and 7e(., .)x for 
fixed scaling parameter X k = 1 are compared, where both give rise to a mesh inde¬ 
pendent method and the latter leads to a large speed up. Note that the choice of 
(., .)x with X k = (ye) -1 leads to the same iterates than choosing ye(., -)x and X k = 1. 
Furthermore, it is discussed in |3j that the choice of je(., .)x can be motivated using 
j"( cp) or by the fact that for the minimizers {^> £ } e> o the Ginzburg-Landau energy 
converges to the perimeter as £ -> 0 and hence yeiy’ellx ~ const independent of 
e « 1. However, since this holds only for the iterates <p k when the phases are sepa¬ 
rated and the interfaces are present with thickness proportional to e, we suggest to 
adopt X k in accordance to this. As updating strategy for X k the following method 
is applied: Start with Ao = 0.005(7e) -1 , then if a k -1 = 1 set X k = X k -i/0.75, else 
X k = 0.75Afc_i and X k = max{A TO j n , min{A mo3; , A*,}}. The last adjustment yields that 
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(A12) is fulfilled. Numerical experiments in [3] show that this in fact produces for 
the choice (., .)x a scaling with A& « (7s) -1 for large k. 

In (3] [5] the effect of obtaining various local minima of the nonconvex optimization 
problem (13)-(15) by choosing different initial guesses (p 0 can be seen. However also 
the other parameters have an influence. 

In this paper we concentrate on comparing different choices of the inner products 
a*; and use herefor the cantilever beam described in [3] with -0o(¥>) = |(1 -if-ip) 
and a quadratic interpolation of the stiffness tensors C(ip). The computation are 
performed on a personal computer with 3GHz and 4GB RAM. First we discuss 
the choice of (,,.) L 2 versus (.,.)x- The choice of the L 2 -inner product leads to 
the commonly used projected L 2 -gradient method. However, (., ,) L 2 does not fulfill 
the assumptions of the VMPT method, since j is not differentiable in L 2 (Q) N or 
L 2 (Q) N n L°°(kl) N . Thus, global convergence is given for the discretized, finite di¬ 
mensional problem but not in the continuous setting. This leads in contrast to the 
choice of (,,.)x to mesh dependent iteration numbers for the L 2 -gradient method, 
which can be seen in Table [TJ The values in Table [T| were computed for different 
uniform mesh sizes h with the parameters e = 0.04, 7 = 0.5, tp 0 = m and tol = 10 -5 
for the stopping criterion V^Il 2 ^ tol. The behaviour of iteration numbers is 
in accordance to our analytical results in function spaces considering h -+ 0. Fur¬ 
thermore, numerical results not listed here show that we obtain for (., ,)x and large 
k scalings A*, « (7s) -1 independent of the mesh parameter h, whereas the L 2 -irmer 
product produces A*, scaled with h 2 . Since the algorithm using the L 2 -inner prod¬ 
uct is equivalent to the explicit time discretization of the L 2 -gradient flow, i.e. of 
the Allen-Cahn variational inequality coupled with elasticity, with time step size 
At = Afc, the scaling A& = 0(h 2 ) reflects the known stability condition At = 0(h 2 ) 
for explicit time discretizations of parabolic equations. 


h 

2-4 

2" 5 

2-6 

2- 7 

2-8 

(-,-)* 

323 

5015 

18200 

57630 

172621 

(■> -)x 

111 

407 

320 

275 

269 


Table 1: Comparison of iteration numbers for (.,.)l 2 and (.,.)x- 


Next we compare (.,.)x with a*, given in (21), which incorporates second order 
information. As experiment we again use the cantilever beam in [3], now with 
e = 0.001, 7 = 0.002, tol = 10 -4 and random initial guess ip Q together with an adaptive 
mesh, which is fine on the interface with h ma x = 2 -6 and h m i n = 2~ n . The parameter 
Afc is updated as described above. The computational costs of one iteration with 
ak given in (21) is significantly higher, since the calculation of Vk(¥>k) requires the 
solution of a quadratic optimization problem with cp e and in addition with 
the linearized state equation (20) as constraints. However, in each PDAS iteration 
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Figure 1: Local minima for the cantilever beam. 


solving the subproblem for fixed k, only the right hand side of (20) changes, namely 
only p. We factorize the matrix in the discrete equation once such that for each 
p only a cheap forward and backward substitution has to be done. In Table [2] the 
corresponding iteration numbers, the total CPU time, the values of the combined 
cost functional as well as of the parts, i.e. the mean compliance and the 

Ginzburg-Landau energy are listed. One observes the drastic reduction in iteration 
numbers using second order information. Due to the mentioned higher costs of 
calculating the search directions the total CPU-time is only halved. Nevertheless, 
this can be possibly improved using a more sophisticated solver for Vki^Pk)- ^ can 
be also observed that the cost and the probably more interesting value of the 

mean compliance is lower. Hence, the different inner products result in different 
local minima, which are shown in Figure [TJ The inner product given in (21) yields 
a finer structure. Also in other experiments we observed a local minima with lower 
cost value for this choice of a*,. 


inner product 

iterations 

CPU time 


fr a 9 ' u* 

EM 

(■> -)x 
a k in (21 

) 

11189 

851 

42h 12min 

19h 

15.07 

14.99 

15.03 

14.93 

20.79 

30.12 


Table 2: Comparison of two different inner products. 


We successfully applied also an L-BFGS update in function spaces (see e.g. 
for the unconstrained case in Hilbert space) of the metric a k , i.e. starting 
a 0 (u,v ) = 7 e(u,v)x we use the update 


a k+1 (u,v) = a k (u,v) - 


a fc(Pfc> M W(Pfc>U + {yk,u),{y k ,v) 
a k(p k ’Pk) (Vk,Pk) 


ra 

with 


in case that (y k ,p k ) > 0, where p k := < p k+1 - ip k and y k := j'(ip k+x )~j'(<p k ), which 
performs very good especially for small 7. Note that - as in the finite dimensional 
case - assumption (A8) cannot be shown for this sequence of inner products, but 
numerical experiments show that the discretized method is mesh independent, see 
Table [3] for the above cantilever beam example, where the maximal recursion depth 
is set to 10. 
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h 

2~ 5 

2-6 

2" 7 

H 1 -BFGS iterations 

85 

88 

86 


Table 3: Mesh independent iteration numbers for the H 1 -BFGS method. 


The following compliant mechanism problem 


1 

mm - 
2 


f (l~ip N )\u-u n \ 2 + 


where the elasticity equation (14) and the constraints (15) have to hold, is more 
difficult. In our numerical analysis the solution process is more sensitive to the 
choice of a*,. Here the above H l -BFGS approach enables us to solve the problem 
in an acceptable time. Until 7£|Vnfc|i,2 < tol = 10 -4 the calculation of the material 
distribution in Figure 2a took 22 hours. It aims to crunch a nut in the middle of 
the left boundary when the force acts on the right hand side from above and below 
and the mechanism is supplied on the left boundary. 

Moreover, we also successfully applied the VMPT method on the following drag 
minimization problem of the Stokes flow using a phase field approach, which is 
analysed in [16]: 


r 1 1 

mhl Jn 2 1 l Vw | 2+ 2 w(v)|u | 2 + 

f a e (ip)uv + f Vu • Vv = 0 Vv e Hl dlv (il) 
u|dn = (1,0) T , J~ ( P = 0.75, -l<ip<l. 

We applied a nested approach in h and e as well as an adaptive grid. As inner 
products we used the above //'-BFGS method and obtained the result in Figure [2b] 
with 188 iterations to obtain tol = 10 -3 , which took 17 minutes. 

A different type of optimization problem is the inverse problem for a discontinuous 
diffusion coefficient, where the discontinuous coefficient a is smoothed by a phase 
field approach and no mass conservation is used mt 

min i [ |« - Moisl 2 + jE(tp) 

s.t. / a(<p)Vw-V£= / gi e H 1 and u- u 0 b s , -1 < <p < 1. 

Jn Jr Jn Jn 

We choose u Q b s as solution of the state equation for ip shown in the upper part of 
Figure [2c] with added noise of 5% and obtain the solution shown in the lower part 
of Figure [2c} 

The VMPT method can also be used for image inpainting using a phase field ap¬ 
proach by considering 


min i||(p-/||^ (nxD) + 7 U((p) 
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(a) Crunching mechanism, (b) Obstacle minimizing drag. (c) Identified coefficient. 
Figure 2: Successful applications of the VMPT method. 


such that (p fulfills (15), where / is the given image and the inpainting is performed 
in D [9j. The method can adjust to the chosen metric II (Q \ D ) and for this problem 
a line search with exact step length can be applied 


The last four mentioned application examples are preliminary results and are under 
further studies. To our knowledge the VMPT-method outperforms the existing 
applied optimization algorithms in these cases. 
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